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Abstract 

We  identify  mechanisms  for  particle  transport  across  a  cross-field  sheath.  We  present 
a  study  of  E  x  B  drift  motion  in  a  vortex  in  which  the  ion  drifts  are  perturbed  by  their 
finite  gyroradii  and  electron  drifts  are  perturbed  by  one  or  more  traveling  waves.  Large 
scale  vortices,  which  are  the  result  of  nonlinear  saturation  of  the  Kelvin-Helmholtz 
instability  resulting  from  shear  in  the  E  x  B  drift  velocity,  have  been  observed  in  plasma 
simulations  of  the  cross-field  sheath1-3.  Small  scale  turbulence  is  also  present,  and  ions 
and  electrons  are  transported  across  the  sheath.  A  vortex  alone  does  not  allow  for 
otPresent  address:  Plasma  Physics  Laboratory,  Princeton  University,  Princeton,  New  Jersey  08543. 
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the  observed  electron  transport  because  the  electron  drift  orbits  simply  circulate.  On 
the  other  hand,  the  ion  motion  can  be  stochastic  from  resonant  interaction  between 
harmonics  of  the  the  drift  motion  and  the  gyromotion,  independent  of  the  background 
turbulence.  The  fluctuations  in  the  ion  density  would  then  give  rise  to  small  amplitude 
wave  spectrum.  The  combined  action  of  the  vortex  fields  and  traveling  wave  fields 
on  the  electron  motion  can  then  lead  to  stochastic  electron  diffusion.  We  study  these 
effects,  showing  that  the  values  of  vortex  fields  observed  in  the  simulation  are  sufficient 
to  lead  to  both  ion  and  electron  stochasticity.  Furthermore,  the  rate  of  the  the  resulting 
diffusion  is  sufficient  to  account  for  the  diffusion  observed  in  the  simulation. 


1  Introduction 

In  previous  studies  of  the  cross-field  sheath,  the  Kelvin- Helmholtz  (K-H)  instability  is  found 
to  saturate  into  large  scale  vortices,  and  in  addition  to  the  circular  vortex  flow,  there  is  a 
bigger  amplitude,  small  scale  turbulence1-3.  The  driving  mechanism  for  the  K-H  instability 
is  the  E  x  B  sheared  flow  that  occurs  due  to  nonuniform  electric  fields  in  the  sheath  near  a 
conducting,  absorbing  wall.  With  the  magnetic  field  taken  parallel  to  the  wall,  the  electric 
field  is  due  to  the  weaker  magnetization  of  the  ions  relative  to  the  electrons,  causing  a  net 
positive  wall  charge.  This  electric  field  is  nonuniform  (large  at  the  wall,  becoming  small  in 
the  plasma),  which  causes  a  significant  velocity  shear;  the  mean  velocity  is  r0  ~  where 
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vt,  is  the  ion  thermal  velocity,  and  the  shear  length  is  L,  ~  5 p,,  where  p,  is  the  thermal  ion 
gyroradius.  A  characterization  of  the  physical  system  is  shown  in  Figure  1.  Our  objective  in 
this  paper  is  to  understand  how  the  vortices  observed  in  the  previous  Refs.  1-3  can  provide 
mechanisms  for  transport.  The  calculation  presented  here  may  be  relevant  to  K-H  vortices 
driven  by  other  mechanisms  as  well.4,5 

There  has  been  considerable  research  on  chaotic  single  particle  motion  in  plasmas,  in 
which  chaos  in  phase  space  has  been  generated  by  the  resonant  interaction  of  the  particle 
gyration  with  magnetic  field  spatial  periodicities6,7,  with  periodic  time  varying  fields8,9, 
and  with  waves10,11.  The  resonant  mechanism  for  the  generation  of  the  stochasticity  and 
the  resulting  particle  diffusion  is  the  overlap  of  neighboring  harmonic  resonances  in  the 
action  space  as  reviewed  in  Refs.  12  and  13.  There  has  also  been  interest  in  stochastic 
E  x  B  motion  and  its  role  in  transport  due  to  drift  waves  both  when  the  stochasticity 
is  caused  by  an  interaction  of  one  or  two  waves14-16  or  when  the  wave  spectrum  itself  is 
turbulent17,18,19.  Our  work  is  in  the  same  spirit  as  Refs.  14-16,  but  our  physical  system, 
and  associated  Hamiltonian  are  quite  different. 

We  postulate  that  the  basic  driving  mechanism  for  the  transport  is  the  following.  The 
nonlinear  ExB  motion  within  the  vortices  generate  harmonics  of  the  vortex  frequency 
that  resonate  with  the  ion  gyrofrequency.  These  harmonic  resonances  generate  secondary 
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resonances  whose  interaction  make  the  ion  motion  stochastic.  The  resulting  interception  of 
ion  orbits  by  the  walls  leads  to  macroscopic  charge  fluctuations  that  generate  a  wave  spec¬ 
trum.  The  waves,  in  turn  have  frequencies  that  resonate  with  the  ExB  drift  frequency  of 
the  vortex  motion  of  the  electrons  to  produce  stochasticity  and  therefore  electron  diffusion. 
The  sheath  physics  is  therefore  of  key  importance  in  setting  up  a  nonuniform  electric  field, 
which  leads  to  vortex  flow  and  therefore  a  frequency  spectrum  associated  with  the  vortex 
motion. 

Although  the  driving  terms  for  generating  the  stochasticity  proceed  from  ions  to  elec¬ 
trons,  the  electron  dynamics  are  easier  to  treat,  so  we  reverse  the  order  of  presentation. 
In  Part  I,  we  concentrate  on  the  electron  E  X  B  motion,  which  is  also  applicable  to  the 
zero  order  ion  drift  motion.  In  Section  2,  we  develop  the  interaction  equations  between  an 
applied  vortex  field  and  one  or  more  wave  fields.  In  Section  3,  we  study  orbits  obtained 
from  numerical  integration  of  the  equations  of  motion.  We  show  results  using  single  and 
multiple  waves,  and  compare  with  the  analytic  predictions.  We  also  show  the  importance 
of  the  edge  velocity  shear.  In  Section  4,  using  a  resonance  overlap  criterion,  we  calculate 
how  large  the  perturbing  wave  needs  to  be  to  cause  large  scale  stochasticity. 

In  Part  II  we  consider  the  interaction  of  the  finite  ion  gyromotion  with  the  vortex.  In 
Section  5,  we  introduce  the  interaction  Hamiltonian  and  the  appropriate  transformations. 
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The  numerical  interactions  and  comparison  with  theory  are  presented  in  Section  6.  In 


Section  7  we  again  use  the  overlap  criterion  to  determine  the  parameter  region  for  large 
scale  stochasticity.  In  Section  8  we  calculate  and  measure  the  diffusion  coefficient  arising 
from  the  ion  stochastic  motion,  and  compare  it  to  a  simulation  of  the  ion  motion  in  a  applied 
vortex  field. 

We  emphasize  that  these  calculations  are  not  self-consistent  as  both  the  vortex  field 
and  the  wave  field  are  imposed.  However,  the  amplitudes  of  those  fields  are  taken  from 
the  self-consistent  simulations,  and  the  results  of  the  calculations  show  electron  and  ion 
diffusion  consistent  with  the  observed  self-consistent  diffusion.  Thus,  we  believe  that  these 
calculations  uncover  the  mechanisms  operating  to  produce  the  self-consistent  diffusion. 

Part  I 

2  Electron  Drift  Motion  Interacting  with  Traveling  Waves 

In  all  calculations  we  will  use  the  configuration  specified  in  Refs.  1-3  which  is  given  in  Figure 
1.  The  magnetic  field  B  is  constant  and  perpendicular  to  the  x-y  plane  (parallel  to  the  wall), 
and  the  electric  field  E  is  in  the  x-y  plane.  We  assume  that  the  time  dependent  variation  of 
the  electric  field  is  much  smaller  than  the  electron  cyclotron  frequency  (u;  <  <*>,*).  In  Part  I 
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we  also  assume  that  all  the  field  spatial  scale  lengths  are  much  bigger  than  the  electron 
gyroradius  Cl).  Hence,  in  this  section  we  can  use  the  zero  order  drift  equation  for 

the  motion  of  the  electrons 


v  = 


E  x  B 

~W~' 


(1) 


This  is  also  applicable  to  the  zero  order  ion  drift  motion.  However,  as  will  be  discussed  in 
Part  II,  finite  gyroradius  and  gyrofrequency  effects  are  also  important  for  ion  motion  in  the 


vortex.  Written  out  in  component  form,  in  terms  of  the  potential  defined  by  E  =  —  V$,  we 


see  that  Eq.  (1)  is  in  Hamiltonian  form  with  x  being  the  canonical  momentum  and  y  the 
position, 


I<9£ 

Bdy' 


V  =  15 


l££ 

B  dx  ' 


(2) 

(3) 


where  $  plays  the  part  of  the  system  Hamiltonian.  For  $  independent  of  time,  the  particles 


follow  the  equipotential  contours  $(x,y)  =  §  [x(f  =  0),y(f  =  0)]. 


We  now  study  the  particle  motion  in  a  time  dependent  electrostatic  potential  given  by 


$(x,  y,t)  =  -B  vq  x  - 


2  L, 


-•»[(=)’♦  (HK 


(4) 


+  e$o]T'4"c,(fcnV~u,nt)> 


where  u0  the  drift  velocity  of  the  vortex  and  L,  is  the  velocity  shear  length.  The  potential 
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has  a  magnitude  $o  and  a  “shape  function”  rj)  specifying  the  shape  of  the  vortex  equipo- 
tentials  that  specifies  the  vortex  flow;  a  and  (3  specify  the  elongation  in  the  x,  y  directions. 
Defining  7  =  then  7  and  e  are  dimensionless  parameters  specifying  the  strengths  of 
the  vortex  and  waves,  respectively.  The  elliptical  shape  for  the  zero  order  contours  in  Eq. 
(4)  was  assumed:  (1)  because  of  reasonable  comparison  with  results2,3;  (2)  for  clarity  and 
ease  of  the  calculations.  The  results  presented  below  can  be  used  for  other  closed  contour 
shapes.  A  contour  plot  of  Eq.  (4)  is  shown  in  Figure  2  with  Gaussian  rp(u)  =  e~7u,  with 
5=1,  and  the  dimensionless  parameters  the  same  as  used  in  the  simulation2,3  7  =  2.5, 
a  =  1.25,  /3  =  3.83,  uo  =  0.44uTi,  p,/I„  =  0.8,  and  kpi  =  3.83.  We  note  that  in  the 
simulations  p,  and  vt,  are  routinely  set  equal  to  one,  such  that  the  dimensionless  quantities 
here  correspond  to  the  quantities  in  the  simulation. 

The  notation  can  be  simplified  by  making  the  following  coordinate  transformation  to 
dimensionless  coordinates  in  a  frame  moving  with  the  vortex 


x 

y 


x 

ap,' 

y  +  v0t 

0P,  ' 


where  fl  is  the  thermal  ion  gyrofrequency.  Dropping  the  overbars  for  notational  convenience, 


the  equations  of  motion  are  now 


X  =  — S-, 

dy 

(5) 

y  =  aT 

(6) 

with  the  transformed  Hamiltonian 

$(x,y,t)  =  -i>  ( 

X2  +  y2)  +  6x2  +  e  y^Ane^-^l 

n 

(7) 

|| 

where  6  =  - 2-J-,  parameterizes  the  amount  of  velocity  shear.  We  will 

2')VTiLa 

assume  6  small; 

using  the  early  values  produces  6 

~  0.1.  The  transformed  k,  Q  in  terms  of  the  laboratory 

frame  k,  u,  are 

kn 

—  0knp,, 

(8) 

=  — ft-1(u>n  +  knV0). 
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Next,  we  make  the  transformation  from  (i,  y)  coordinates  to  their  corresponding  action- 
angle  coordinates  (J,9),  using  the  partial  generating  function 


F\{yJo)  =  -y2  cot0o. 


The  new  and  old  variables  are  related  by 


x 


Jo 


dF\ 

dy  ’ 

dFx 

d8o' 


(10) 


(11) 


8 


(12) 


which  after  rearranging  gives 


y/2JoCOsdo, 

(13) 

\/2Jq  sin  9q. 

(14) 

For  drift  motion  in  a  quadratic  well,  J o  gives  the  area,  divided  by  2x,  inside  the  curves  of 
constant  Hamiltonian,  and  do  is  a  uniformly  rotating  angle.  Defining  the  unchanged  wave 
phases 


9n  =  unt,  n  =  1,2,3, 


(15) 


It  is  necessary  to  increase  the  system’s  degrees  of  freedom  to  keep  the  angle  variables 
periodic.  We  have  assumed  the  waves  are  independent  (no  harmonics).  We  now  obtain  the 
Hamiltonian  H  —  E o  +  (E\ , 

n=N 

Ho(Jq,Ji,J2,—,Jn)  =  -1P(2Jo)+  ^  &nJn,  (16) 

n=  1 

H\(Jo,9o,9udt,  —  -,9n)  =  -2Jocos2^o  (17) 

Af  l—+oo 

+  E  JiBCknV2J o)e,(,e°-H 

n=  1  /=— oo 

where  JtB  are  Bessel  functions  of  integer  order,  which  arise  from  the  expansion  of  (14)  in 
the  exponent12,13.  From  this  form,  it  is  seen  that  each  wave  adds  one  degree  of  freedom  to 
the  system,  but  a  single  wave  generates  an  infinite  set  of  resonances  between  harmonics  of 
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the  drift  frequency  and  harmonics  of  the  wave  frequencies.  The  second  term  in  Eq.  (16) 
preserves  the  canonical  form  for  n  >  1.  Eqs.  (13)-(17)  make  no  reference  to  the  assumed 
shape  of  the  vortex  contours,  beyond  the  form  indicated  in  (7).  Hence,  the  results  here  and 
following  can  be  generalized  to  other  vortex  shapes  satisfying  this  form.  In  a  general  case 
J0  would  be  the  area  inside  the  closed  contour  of  the  potential  divided  by  2tt  (with  e  set  to 
zero)16. 

For  (  equal  to  zero,  Jo  is  a  constant  of  motion  and  the  drift  frequency  around  the  vortex 
Ud  is  defined  by: 


J>d  —  & 0  — 


dip(2J0) 
dJ0  ' 


(18) 


For  small  e,  we  expect  resonances  when  lu>d  =  Tnun,  where  /,m  are  integers.  If  ip  is 
nonlinear,  which  is  the  case  of  interest,  then  £Jd  is  a  function  of  Jo  and  there  will  be  many 
resonances  for  a  given  u„.  These  resonances  generate  islands  in  the  phase  space  which  can 
overlap  to  provide  a  mechanism  for  an  electron  to  stochastically  change  J o,  hence  providing 


a  mechanism  for  transport6,13. 


3  Numerical  Results  for  E  x  B  Motion 

In  the  K-H  simulations,  small  amplitude  and  short  scale  waves  traveling  in  the  y-direction 
were  observed  throughout  the  simulation  (  e  <  1,  k pi  ~  1  and  ~  1  ).  In  this  section  we 
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show  how  the  short  wavelength  and  small  amplitude  waves  can  provide  a  mechanism  for 


transport.  We  add  a  small  amplitude  wave  component  to  the  large  scale  vortex  flow  and 
study  the  resulting  E  x  B  particle  motion.  A  single  wave  of  appropriate  amplitude  should  be 
sufficient  to  obtain  global  stochasticity  and  we  concentrate  on  this  case.  We  again  choose: 
i’(u)  =  —  e“5u,  and  we  neglect  velocity  shear  (6  =  0).  The  potential  is  then 

$(x,y,t)  =  -  exp  ^"(x2  +  I/2))  +  * cos(ky  -  ut)  (19) 

In  this  case,  with  e  =  0,  the  drift  frequency  is:  uj  =  exp(- Jo)-  The  values  for  a,  /?,  and 
7  are  approximately  1.25,  3.83  and  2.5  for  the  simulation  presented  in  Refs.  [2,3].  In  the 
norm-dized  units  the  maximum  drift  frequency  is:  Cjb<max  =  1,  which  in  the  laboratory  units 
is*  id b,max  —  ^b,max  ^  0.5212.  This  compares  reasonably  well  to  the  observed  value  of 
0.512  [3]. 

We  integrate  the  equations  of  motion  Eq.  (5)  and  (6),  for  a  set  of  orbits  using  the 
classical  fourth  order  Runge-Kutta  method20.  The  wave  potential  is  given  by  Eq.  (19)  with 
k  =  3.83,  u  =  0.96.  Our  choice  of  k,u  was  made  so  that  kp,  =  1  and  u>  was  a  representative 
point  from  the  power  spectrum  at  kp,  as  1.  Fig.  8  in  Ref.  3  shows  the  power  spectrum 
peaked  at  u>0  ss  -kv o,  where  vo  ~  0.44uxi,  but  is  fairly  broad  at  kpi  k  1  with  a  range 
between  -2 u>0  and  zero.  In  Figure  3  we  display  integrations  from  a  set  of  initial  locations: 
x i(t  =  0)  =  0.05  t  and  y,{t  =  0)  =  0  with  i  =  1,2,  ...,20.  Surface  of  section  plots  in  (x,  y )  are 
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shown  for  1000  periods  of  the  perturbing  wave  (so  that  there  are  20,000  total  points).  The 
surface  of  section  is  defined  at  constant  phase  of  the  wave  variable,  i.e.,  when  ut  —  2nir , 
where  n  integer.  If  a  constant  of  the  motion  exists,  then  it  has  been  generally  shown13  that 
i  and  y  will  lie  on  a  smooth  curve  in  this  surface.  If  the  constant  does  not  exist,  then  the 
trajectory  can  puncture  the  surface  of  section  over  a  two-dimensional  area. 

We  observe  the  transition  between  these  two  types  of  behavior  in  Figs.  3(a)  and  3(b). 
The  plots  show  the  change  in  character  as  6  is  varied  over  a  range  consistent  with  the  wave 
amplitudes  observed  in  the  simulation,  (a)  e  =  0.015,  and  (b)  e  =  0.02.  In  Fig.  3(a)  the 
motion  is  mainly  confined  to  smooth  curves  or  bounded  in  small  areas  between  smooth 
curves.  The  latter  regions  develop  near  the  separatrices  of  island  chains  formed  by  reso¬ 
nances  between  the  drift  motion  around  the  vortex  and  the  drift  motion  in  the  wave  field. 
In  these  bounded  separatrix  layers  the  constant  of  the  motion  does  not  exist,  but  this  has 
little  effect  on  the  overall  dynamics.  In  Fig.  3(b),  the  separatrix  layers  of  neighboring  islands 
have  overlapped  to  make  a  broad  layer  of  stochasticity.  This  is  clearly  seen  in  the  escape  of 
trajectories  to  much  larger  values  of  the  coordinates.  If  trajectories  intersect  sources  and 
sinks  of  plasma  they  would  then  result  in  diffusion.  We  present  a  quantitative  calculation  of 
diffusion  in  the  next  section.  Figure  4  is  the  surface  of  section  plots  using  the  action-angle 
coordinates  (Jo, 9 o)  for  the  (  —  0.015  case. 
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If  we  consider  that  more  than  one  wave  is  in  the  perturbing  field,  the  potential  would 
have  the  following  form 

1  'i  M 

-~(x2  +  y2)  [  +  «  £  Am  cos (km  -  umt  +  0m),  (20) 

We  would  expect  that,  for  the  same  overall  wave  amplitude  level,  phase  decorrelations 
might  increase  the  level  of  stochasticity.  As  an  example,  we  use  four  waves  ( N  =  4), 
two  with  positive  w,  two  with  negative  w,  all  having  kpi  =  1.024,  and  the  amplitudes 
representative  of  the  power  spectrum  given  in  Ref.  3,  Figure  8.  The  phases  0m  were  chosen 
to  be  random.  The  parameters  used  for  the  four  waves  n  =  (1,2, 3, 4)  are:  k„  =  4.16, 
lin  =  (-0.6, -0.3, 0.4, 0.7),  eAn  =  (0.004,0.006,0.01,0.012),  and  0m  =  (1.1, 5.5, 3.0, 0.0). 
Figure  5  shows  a  surface  of  section  plot,  where  points  are  plotted  when  u ;4f  =  m27T,  m  = 
0,1,2, 3,...  This  surface  of  section  is  chosen  because  mode  four  has  the  largest  amplitude, 
so  that  residual  structure  can  still  be  seen.  The  lack  of  a  true  surface  of  section  in  which  to 
see  regular  motion  prevents  a  quantitative  comparison  with  the  single  wave  case,  but  the 
impression  is  of  increasing  stochasticity  with  more  than  one  wave. 

We  can  also  investigate  the  role  velocity  shear  plays  in  the  E  x  B  particle  motion.  We 
do  this  by  adding  bx2  to  Eq.  (19).  The  velocity  shear  permits  plasma  to  flow  past  and 
interact  with  the  vortex.  This  type  of  shear  is  observed  in  the  self-consistent  simulations. 
The  velocity  is  strongest  at  the  wall  and  drops  off  to  zero  towards  the  center  of  the  plasma. 


$(x,y,t)  =  -  exp  | 
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The  form  x2  used  for  modeling  the  shear  only  applies  up  to  x  =  2 La/(ap,),  at  which  point 
the  drift  velocity  goes  to  zero,  as  can  be  seen  from  Eq.  (4).  The  velocity  shear  and  vortex 
provide  a  “snow  shovel”  mechanism,  by  which  the  vortex  drifts  along  the  wall  and  scrapes 
off  interior  plasma  which  is  within  the  reach  of  the  vortex  (~  5 p,).  This  mechanism  is 
shown  in  Figure  6,  using  the  example  of  c  =  0.03,  an  6  =  0.015  for  a  single  wave.  Four  test 
particle  orbits  are  shown  just  at  the  edge  of  the  vortex. 

4  Island  Formation  and  Overlap  Criterion  for  Global  S- 
tochasticity 

We  begin  by  studying  a  single  wave.  To  find  out  for  what  range  of  parameters  global  (large 
scale)  stochasticity  occurs,  we  calculate  the  overlap  of  the  first  two  integer  resonances: 
tJd  =  w/m,  m= 1,2.  These  two  resonances  are  chosen  because  we  are  interested  in  waves 
with  Cj  *>d,maxi  and  these  are  the  primary  resonances  observed  in  the  simulation  results 
presented  in  Section  3.  We  make  a  transformation  to  the  slowly  varying  phase  close  to 
resonance:  4>0  =  mOo  -  9\,  and  the  fast  phase  of  interest:  =  9i/m.  Using  the  generating 

function: 

^(/o,  It;  Bo,  0i )  =  {m0 o  -  9\  )Iq  4-  Si Ii/m, 
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we  obtain: 


J0  =  J0/m ,  (21) 

Jj  =  Jo  +  mJ\ , 

<£o  =  tn^o  —  $i , 

4>i  =  Oi/m. 

Which  explicitly  exhibits  the  slow  phase  <f>o.  The  new  Hamiltonian  is  then  H  =  Hq  +  cJTi 
where 

Ho(IoJx)  =  -4>(2ml0)  +  w(Ji/m  -  J0),  (22) 

^  ( Jo,  0o,^i )  =  -2mjocos2(<£o/m  +  </>i)  (23) 

f=+ oo 

+  €  Ji8(fcy2^)e’^+('-,7,)^), 

/=  — OO 

and  we  have  taken  Aj  =  1.  Next,  we  average  over  the  fast  phase  4>\  to  obtain: 

H  =  -rp(2mIo)  -  CjIq  +  6mIo  +  f  J®(fcv/,2mJ0)e,,*>0  (24) 


where:  H  =  ^  Hdtfri.  At  this  point  we  drop  the  zero  subscript  since  the  the  phase 


,  .  .  dH  dH 

averaged  system  is  now  one  dimensional.  Setting  — —  =  0,  — -  =  0,  gives  the  fixed  points, 

01  dtp 


which  are 


Ifp  — 

l{-u/2m). 

(25) 

<*7p  1.2  = 

(0,?r), 

(26) 
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where  we  have  neglected  terms  of  order  e  in  determining  7/p.  We  can  then  expand  I  = 


Ifp  +  A 7  to  obtain13: 


A/2 

A 27  =  G—  -  2^  COS(j> 


where  G  =  4m2^”(2m//p),  F  =  —tjB(ky/2TnIjp).  The  maximum  A/  is 


f¥_  r-eJ^ygmTg 

V  G  m2t/>"(2m/yp) 


As  an  estimate  of  the  onset  of  large  scale  stochasticity  we  use  the  “two-thirds  rule”,  which 
corresponds  to  the  destruction  of  the  last  KAM  surface  between  islands(at  K=1  for  the 


standard  map)13: 


A/maJ(m  —  1)  4~  A7mai(m  —  2)  ^  2 
Ifp(m  =  2)  -  IfP(m  =  1)  _  3 


To  proceed  further  the  vortex  shape  needs  to  be  specified.  We  use 


0(u)  =  e 


which  gives  a  qualitative  fit  to  the  vortex  observed  in  Refs.  2  and  3.  Using  this  shape  we 


calculate  the  location  of  the  mth  resonance  in  action  space  to  be  at 


Ifp  -  -  —  ln(w/m). 


Using  (30)  to  calculate  in  (28),  the  approximate  width  of  the  mth  resonance  is 


A7mai  =  2 


(jB(ky/2mIfp) 
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We  can  now  use  Eq.  (29)  with  Eqs.  (31)  and  (32)  to  calculate  the  critical  value  for  the 
perturbation  strength: 

N-2171^  +\]^ ■  (33) 

The  value  of  ec  gives  an  estimate  of  the  strength  of  the  background  perturbing  wave 
which  is  needed  to  cause  large  scale  stochasticity.  The  predicted  critical  value  for  the  per¬ 
turbation  strength  obtained  using  Eq.  (33)  is:  ec  =  0.012,  which  is  in  reasonable  agreement 
with  the  numerical  results  which  show  the  transition  to  large  stochasticity  occurring  be¬ 
tween  Figs.  3(a)  and  3(b).  If  we  assume  that  « is  sufficiently  large  that  the  resonance  islands 
strongly  overlap,  i.e.  (29)  is  well  satisfied,  there  is  good  mixing  of  the  particles  in  the  vortex 
and  we  can  then  estimate  the  stochastic  diffusion  coefficient  by  its  quasilinear  value 

D.  ~  AL^  (34) 

where  Qp  is  the  drift  frequency  and  AL  a  characteristic  diffusion  distance  in  physical 
space  expressed  in  dimensional  units.  We  associate  AL  with  the  half  island  size,  given  in 
(32),  where  AL  and  AImax  are  related  through  the  transformation  equations.  Substituting 
Eqs.  (32),  (31)  and  (18)  into  (34)  we  obtain,  in  dimensional  units 


We  note  that  the  factor  =  pffl  restore  the  dimensionality,  but  does  not  imply  that  D, 
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has  Bohm  scaling,  since  the  dimensionless  factors  can  also  change  with  this  ratio.  We  will 
compare  this  result  with  the  ion  diffusion  obtained  in  Sec  II.  From  Eq.  (35)  we  obtain  the 
value  of  the  stochastic  diffusion  coefficient  is  Da  ~  0.07^-.  It  is  worth  pointing  out  that 
the  similar  issue  has  been  addressed  in  a  quite  different  system  in  which  equilibrium  and 
perturbed  states  are  two  dimensional  periodic  flows  with  perturbed  flow  propagating  in 
y-direction.21 

Part  II 

5  Interaction  of  the  Ion  Drift  and  Gyromotion 

Because  the  gyro-orbits  of  the  ions  are  comparable  in  size  to  the  vortex,  the  drift  approx¬ 
imation  can  no  longer  be  used.  Furthermore,  the  gyromotion  adds  a  second  degree  of 
freedom  which  allows  resonances  even  in  the  absence  of  a  wave  field.  We  now  consider  the 
formulation  of  the  ion  dynamics  including  the  full  gyro  and  drift  motion. 

In  the  uniform  magnetic  field  B0i  from  the  vector  potential  A0(y)  =  -Boyx,  the  total 
Hamiltonian  is 

B  =  Y\{\P:  +  eBoy\2  +  2 JfPv  +  (36) 
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with  the  electrostatic  potential  well 

${x,y,t)  =  -Bqv0x  -  $0 ip  +ff  (^~“)]  •  (37) 

Here,  ro  is  the  drift  velocity  of  the  vortex,  ip,  f  and  g  are  functions  specifying  the  shape  of  the 
vortex.  For  analysis  of  ion  motion,  we  concentrate  on  the  motion  inside  a  vortex  potential 
well  and  therefore  can  neglect  the  velocity  shear  and  the  perturbed  wave  background  in 
Eq.  (37)  in  contrast  to  Eq.  (4). 

As  in  Part  I,  we  eliminate  t  by  transforming  to  the  frame  moving  with  the  vortex  velocity 
v0:  (x,y,px,py)  -*  (u,v ;PU,PV),  with  new  coordinates  in  dimensionless  units: 


u 


v 

Pu 


X 

■> 

Pi 

y  -f  vpt 
Pt 

M  vj,  vji  ' 


Pv 

H 


Py  Vo_ 

Mv  T,  t'Ti  ’ 

H 

t: 


Following  the  treatment  of  Smith  and  Kaufman,11  we  transform  to  guiding  center  coordi¬ 
nates  (X,Y)  and  local  polar  gyrocoordinates  (p,<p),  (u,  v;Pu,  Pv)  — ►  (Y,X\(p,P^),  by 


Y  =  v  +  psm<p,  (38) 

X  =  u  -  p  cos  <p, 
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<t>  =  tan  *( — - — ), 

*V 

P, ,  =  ip2  =  j[(P.  +  »)2  +  Pi)- 

This  transformation  allows  separation  into  fast  gyromotion  and  slow  drift  motion  so  that 
we  can  use  the  method  of  averaging.  The  Hamiltonian  is  then 

H  =  P9-  7 $[f(X  +  p cos <t>)  +  g(Y  -  psin<t>)}.  (39) 

As  in  Part  I,  we  exhibit  this  structure  by  transforming  to  action- angle-like  variables.  We 
already  have  the  pair  (P#,  4>)  for  the  gyromotion.  The  transformation  of  the  drift  variables 
to  (J,9)  follow  as  in  (10).  Choosing 

X  =  VtJ  cos  9,  (40) 

Y  =  y/2J  sin  9. 

which  is  equivalent  to  Eqs.  (13,14).  The  Hamiltonian  becomes: 

H  =  P4-  7i’[/(v/2 7 cos  9  +  pcos<t>)  +  g(V2Js'm9  -  psin0)].  (41) 


6  Numerical  Results  for  Ion  Motion 


In  order  to  compare  with  2d  simulation  results,  we  choose  two  forms  of  the  potential  well: 
(a),  as  in  Part  I, 


7  0(x,y)=  7  exp 


(42) 
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which  are  integrated  by  a  Boris  mover.22  Since  the  Hamiltonian  is  independent  of  time, 
its  numerical  value  has  been  checked  and  is  essentially  conserved  during  the  integration. 
In  order  that  a  surface  of  section  exhibit  nonintersecting  orbits,  we  start  with  all  particles 
having  the  same  Hamiltonian  and  vary  the  initial  particle  positions  and  velocities.  In 
Figs.  8  and  9  we  display  integration  from  a  set  of  initial  locations  on  the  circle:  x  ,  = 
r0a  cos  (fo, ,  y,  =  r0(3  sin  4>o, ,  vxt  =  uqo  cos  0Qi ,  and  Vy,  -  vqo  sin  #o«  with  <t>ox  and  90 ,  uniformly 
distributed  between  0  and  2ir. 

Because  of  the  complexity  of  the  orbits,  we  visualize  the  motion  both  by  plotting  a 


particle  trajectory  in  physical  space  and  a  number  of  trajectories  in  the  surface  of  section. 
From  the  particle  trajectory,  we  find  the  ratio  of  gyrofrequency  fic  (where  fic  is  the  in¬ 
dividual  particle  gyrofrequency  with  a  specific  velocity)  to  the  E  x  B  drift  frequency  flj 
and  obtain  a  qualitative  feeling  for  the  dynamics.  As  in  Part  I,  to  determine  if  an  orbit  is 
chaotic  or  regular,  we  generate  a  two-dimensional  surface  of  section,  plotting  guiding  center 
coordinates  Y  vs  X  at  a  particular  value  of  <t>  (<t>  =  0  is  convenient)  .  This  is  accomplished 
numerically  by  solving  the  equations  of  motion  in  Cartesian  coordinates  and  calculating  vx 
at  each  time  step;  when  vx  changes  sign  (which  is  equivalent  to  4>  —  0),  we  calculate  X  and 
V'  according  to  Eqs.  (38). 

We  first  demonstrate  the  effect  of  a  resonance  between  a  harmonic  of  the  drift  motion 
and  gyromotion  for  potential  form  (a),  illustrating  the  “divided  phase  space”  of  regular 
and  stochastic  orbits.  Since  the  inflection  point  of  the  vortex  has  zero  velocity  shear,  most 
of  the  particles  in  that  neighborhood  have  small  shear  and  therefore  the  resulting  islands 
have  large  amplitude.  We  therefore  take  the  inflection  point  r0  =  1.0  as  the  initial  position 
and  vqo  =  0.71rre  Fig.  8(a)  shows  a  primary  4:1  resonance  between  the  gyrofrequency 
and  drift  frequency  of  a  particle  orbit  in  configuration  space,  clearly  showing  the  resonance 
of  four  gyro-orbits  to  one  drift  orbit.  Fig.  8(b)  shows  successive  intersections  of  a  number 
of  trajectories  with  the  surface  of  the  section,  with  the  orbit  of  Fig  8(a)  being  one  of 
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the  regular  (island)  orbits.  The  upper  boundary  of  the  phase  space  is  limited  by  energy’ 


conservation.  The  inner  space  is  bounded  by  KAM  curves. 

The  generic  resonance  features  exist  for  various  potential  forms.  Figs.  9(a)-(b)  show  the 
metamorphoses  of  the  Y  -  X  surface  of  section  with  varying  initial  position  tq  and  velocity 
uoo  of  the  potential  form  (b).  For  fixed  initial  position  r o  =  0.5,  Fig.  9(a)  shows  that,  for 
voo  ~  1.55i'7,,  the  motion  is  regular  with  a  5-1  resonance  inside  the  potential  well  and  higher 
resonances  near  the  separatrix.  As  we  allow  t>oo  to  increase  to  1.70ux«i  the  islands  become 
larger  with  the  higher  order  islands  overlapping  to  become  a  seed  for  stochastic  orbits, 
developing  into  a  stochastic  layer  surrounding  the  separatrix,  as  depicted  in  Fig.  9(b).  If 
coo  is  increased  further,  vno  >  1.95 vj,,  some  orbits  can  wander  out  of  the  vortex  cell  and 
escape.  If  the  initial  particle  velocities  are  high  enough,  then  the  motion  near  the  center  of 
the  vortex  can  also  be  stochastic,  as  the  particle  gyroradii  are  sufficiently  large  to  extend 
into  the  nonlinear  part  of  the  vortex  potential.  For  the  parameters  chosen  for  potential  (b). 
particle  motion  can  be  stochastic  and  escape  if  uqq  >  2.72t>x,  even  with  r0  =  0  initially. 
From  Figs.  9(a)-(b),  we  conclude  that  the  motion  becomes  stochastic,  spreading  out  from 
the  separatrix  as  the  velocity  increases. 
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7  Formation  of  Islands  and  an  Overlap  Criterion 


If  the  gyromotion  and  drift  motion  cannot  be  separated  by  averaging,  resonances  that  are 
present  in  the  Hamiltonian  between  gyromotion  and  harmonics  of  the  drift  motion  become 
important.  These  resonances  lead  to  islands  in  the  phase  space  that  have  their  own  local 
phase  space  structure.  Depending  on  the  dynamics,  these  resonances  may  be  well  separated 
or  close  together.  If  closed  together,  then  their  overlap  leads  to  bands  of  stochasticity  result¬ 
ing  in  diffusive  motion  across  the  vortex.  If  well  separated,  the  stochasticity  may  develop 
by  the  interaction  of  second-order  resonances.  These  resonances  are  between  harmonics  of 
the  primary  island  oscillation  and  the  fundamental  drift  frequency  f The  analytical  treat¬ 
ment  of  either  type  of  resonance  can  be  performed  by  the  transformation  to  the  resonant 
frame,  as  done  in  Part  I.  The  new  coordinates  measure  the  slow  oscillation  of  the  variables 
about  their  values  at  resonance,  which  is  an  elliptic  fixed  point  of  the  new  phase  plane. 

In  the  following,  we  proceed  as  in  previous  section  to  obtain  the  motion  near  an  elliptic 
singular  point.  We  choose  the  potential  well  as 


such  that  the  cosine  term  has  a  simple  analytic  expansion  to  exhibit  the  infinite  set  of 
resonances.  We  transform  to  the  action  angle  coordinates  as  in  Eq.  (40)  and  use  the  Bessel 
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function  identity 


OO 

exp(iasin<£)  =  ^  JB(a)exp(tZ</>),  (46) 

/=  —  OO 

to  write  the  Hamiltonian  in  a  form  where  the  resonances  are  explicitly  exhibited: 

E  =  P+  -  -^{(V2Jcos0  +  pcos(f>)2  (47) 

a 1 

—  4j-(\/2Jcos0  +  pcos<j>)4  +  — }  +  ^  ^  {^if  (2v//2J)^(S(--^) 

®  r»=-oo/=— oo 

+  J*(-2VT7)jf(Y)}e,(n<,+/*). 

Noting  the  property  of  the  Bessel  function  JB(-z)  =  (-l)iJjB(i)  and  JB  (x)  =  (-1  ),J/B(i) 
,  we  can  conclude  that  the  resonances  only  have  significant  amplitude  when  n  and  Z  are  both 
odd  or  both  even.  The  strongest  primary  resource  has  the  lowest  order  Bessel  function 
for  which  =  mflj.  Keeping  only  Z  =  ±1  and  n  =  ±m  terms  (where  m  is  odd)  the 
Hamiltonian  which  describes  the  resonance  i-  ^Iven  by 

Hr  =  P#  -  -^{(\/2Jcos0  4- pcos<£)2  (48) 

a 1 

r  _  2  _ 

-  4(\/27 cos  6  +  pcos<t> )4  +  -7-}  -  ^J®(— \Z2J)Jf(— )  cos {m9  -  <f>). 

a*  2  2  00 

Transforming  to  a  slow  variable  near  resonance  by  the  generating  function 

F2  =  (m9  -  <j>)j  +  <j>P<t>,  (49) 

we  have  the  new  variables  as  in  (21)  for  which  the  slow  variable  9  =  m9  -  <b  is  explicitly 
exhibited. 
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Averaging  over  the  rapidly  varying  angle  <£,  the  Hamiltonian  becomes 


S  =  +  (50) 

-  pm(^V/2mi)JiB(y)co^ 

where  Cl  =  1  —  is  the  normalized  gyrofrequency  modified  by  the  vortex  potential,  and  we 
have  assumed  the  distance  from  the  center  of  the  vortex  to  the  fixed  point  R  =  \/2mJ  1- 

The  location  of  the  fixed  points  6fp,  Jfp  are  determined  as  in  (25)  and  (26)  at  §fp  =  0,  tt, 
and 

2mJ» = |k[m  - 1  +  T1-  (51) 

The  island  Hamiltonian  about  the  fixed  points  J )p  is 

AH  =  G(Jfp){-^-~  F(Jfp)  cos6  (52) 


where 


G(J,P) 


d2H 

dP 


%£(*  m)2 


F(Jfp)  =  lj%^yj2mJfp)J?(  ^) 


(53) 


The  island  frequency  near  the  stable  point  is  now 

n,  =  „)  =  (54) 
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and  the  half  island  amplitude  is 


where  the  drift  frequency  near  the  center  of  the  vortex  is  given  by  Eq.  (18)  in  the  dimen¬ 
sionless  units  with  a  =  a,  and  0  =  2b/ir. 


fi,  *  (56) 

We  can  now  calculate  the  ratio  of  the  sum  of  the  two  neighboring  island  widths  and  then 
apply  the  overlap  criterion  as  in  (29).  An  alternative  procedure,  which  has  been  shown  to 
be  equivalent,13  is  to  calculate  the  island  frequency,  via  (54)  and  then  take  the  ratio  to  the 
next  higher  frequency,  which  in  this  case  is  flj.  The  transition  to  large  scale  stochasticity 
in  the  neighborhood  of  the  island  is  then  given  by 


>  1 
nd  ~  i' 


(57) 


Choosing  parameters  a  =  1.59,  b  =  7.5,  m  =  7  and  taking  a  drift  orbit  near  the  y-sepaxatrix 
R,  =  =  7  from  Fig.  9  (calculated  iZ,  value  from  Eq.  (51)  is  smaller  a s  R,  — 

\j2mj fp  ss  3.15),  we  obtain  =  0.12.  This  value,  in  fact,  underestimates  the  ratio,  as  the 
local  value  of  drops  rapidly  when  approaching  the  separatrix.  A  more  exact  calculation 
near  the  separatrix  can  be  performed,  but  requires  more  mathematical  effort.13  We  here 
simply  note  that  the  ratio  obtained  near  the  separatrix  would  satisfy  (57)  if  the  local  value 
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of  ild  were  used.  We  can  do  the  same  calculation  near  the  center  of  the  vortex  for  m  =  5, 
and  Rc  =  \j2rnJjp  =  2  taken  from  Fig.  9,  to  get  =  0.02  The  message  here  is  that  it  is 
easy  to  obtain  stochasticity  near  the  y-separatrix  of  the  vortex  potential  structure  but  the 
orbits  are  mainly  regular  near  the  center  of  the  vortex.  Comparing  these  calculations  with 
the  numerical  results  shown  in  Fig.  9,  we  see  that  the  five-island  chain  has  little  surrounding 
stochasticity,  as  expected.  A  significant  factor  left  out  in  the  drift  frequency  and  the 
fixed  point  JfP  calculations  is  finite  gyToradius  effects,  which  should  give  a  critical  value  of 
velocity  (  or  gyroradius  )  for  the  onset  of  stochasticity,  as  we  see  in  Fig.  9. 

8  Transport  due  to  Ion  Stochastic  Motion 

The  calculation  here  is  motivated  by  self-consistent  simulation  results  from  2-d  bounded 
magnetized  codes  which  have  indicated  that  there  exists  continuous  particle  transport  in  a 
cross-field  plasma  sheath  with  the  large  scale  vortex  potential  structure.  From  dimensional 
considerations,  an  estimate  of  ion  stochastic  diffusion  may  be  made.  During  a  half  gyromo- 
tion  time  xfi-1,  the  particle  is  displaced  over  a  distance  2 p<,  but  the  motion  is  correlated 
such  that  after  a  full  gyroperiod  it  returns  almost  to  its  initial  position  but  with  a  small 
displacement.  If  as  a  result  of  the  ion  gyromotion  resonance  with  E  x  B  drift  motion,  the 
ion  motion  becomes  stochastic  on  this  fast  time  scale,  then  the  successive  displacements  are 
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independent.  A  estimate  of  the  decorrelation  scale  length  is  an  island  size.  With  this  scale 
length,  which  is  that  we  used  in  Part  I,  the  diffusion  is 


(58) 


where  Lj  scales  from  the  island  action  as 


L )  =  2mAJmax 


(59) 


Substituting  Eqs.  (55)  into  (59),  and  (59)  and  (56)  into  Eq.  (58),  we  obtain 

<6o) 

where  we  restore  the  dimensionality  through  the  factor  Slpf  -  jb-  We  compare  this  estimate 
with  the  numerically  determined  diffusion  below.  We  emphasize  here  that  it  is  hard  to 
determine  whether  the  scaling  of  D,  is  Bohm-likein  Eq.  (60)  since  the  parameters  a,  6, 7,f 
and  p  probably  also  change  when  the  factor  f Ip?  =  changes.  We  have  extracted  the 
dimensional  units  for  convenient  comparison  with  the  previous  simulation  studies  by 
Theilhaber  and  Birdsall.3 

In  order  to  numerically  calculate  the  convective  diffusion  coefficient,  it  is  convenient  to 
modify  the  potential  well  with  two  vortex  structures  as 


*  =  (oi) 

apt  ap,  api  api  2  bp, 
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where  H  is  the  step  function.  The  system  is  considered  confined  to  a  domain  -xq  <x  <  xq 
and  —  y0  <  y  <  yo,  where  xq  =  2.5  and  y0  =  b.  Beyond  this  domain,  values  in  y  are 
taken  to  be  periodic.  The  particles  leaving  the  boundaries  at  x  =  ±x<j  are  assumed  to 
lost  to  the  wall  (hence,  there  is  a  sink)  and  reintroduced  at  x=0,  Maxwellian  distributed 
in  velocity  and  y  determined  by  energy  conservation  (  hence,  there  is  a  source  at  x=0). 
When  the  system  reaches  a  steady  state,  the  strength  of  the  source  is  obtained  by  counting 
the  particles  passing  through  x  =  ±x0 .  The  particles  are  initially  randomly  distributed  in 
the  range  x  =  0  and  -yo  <  y  <  yo  and  have  a  Maxwellian  distribution  in  velocity  with 
temperature  T,.  Similar  methods  were  previously  used  in  calculation  of  effective  diffusion 
in  laminar  convective  flows23  and  for  calculating  the  diffusion  in  standard  mapping24.  By 
introducing  a  source  at  x=0,  the  effective  diffusion  coefficient  is  given  by 

D’p  =  5/V2n  ~  5(2x0)2/fV,  (62) 

where  N  is  the  total  number  of  particles  in  the  simulation  region  x  >  0  and  5  is  the  total 
number  of  particles  entering  at  x=0  per  unit  time  from  boundary  x  =  xo- 

In  the  numerical  calculation,  the  system  evolves  for  a  few  gyro-periods  to  reach  a  steady 
state.  Without  the  vortex  potential  well,  the  particles  are  confined  by  magnetic  field,  and 
no  particles  diffuse  to  the  wall.  As  the  amplitude  of  potential  increases,  the  loss  of  the 
particles  increases.  Figure  10(a)  shows  a  plot  of  the  number  of  particles  reaching  the  wall 
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as  a  function  of  time  for  7  =  2.5.  In  this  plot,  2500  particles  were  advanced  for  500  gyro- 
periods.  The  total  number  of  particles  n(t)  escaping  is  plotted  against  t,  and  5  is  found 
to  be  In  this  case,  the  numerical  calculation  of  single  particle  trajectories  found 

D‘p  =  0.03^,  which  is  in  reasonable  agreement  with  the  theoretical  estimate  given  by 
Eq.  (60)  D1/1  =  0.023^.  In  the  theoretical  calculation  we  have  chosen  m  =  7  and  taken  a 
drift  orbit  near  the  y-separatrix  Rt  =  \j2mjjp  =  7  from  Fig.  9.  Comparing  these  results 
to  Theilhaber  and  Birdsall,3  they  obtain  D‘,m  =  0.04^,  is  also  in  reasonable  agreement 
with  the  results  here.  Again,  we  point  out  that  the  simulation  used  pi  =  1  and  vri  =  1,  so 
that  the  ^  scaling  was  not  checked.  Simulations  with  varying  parameters  could  produce 
variations  of  ^  and  the  size  of  the  vortex  in  Eq.  (60). 

In  comparing  the  properties  of  the  surface  of  section  plots  given  in  Section  6  with  the 
diffusion  results,  it  is  necessary  to  understand  the  qualitative  features  of  the  diffusive  motion. 
In  the  surface  of  section  plots,  some  particles  lie  on  regular  orbits  inside  the  vortex  well 
and  other  particles  diffuse  across  the  separatrix.  In  that  case,  we  might  expect  that  the 
density  would  be  uniform  in  the  interior  of  the  vortex,  because  of  the  rapid  drift  motion, 
so  that  the  global  density  would  appear  as  a  steep  gradient  confined  to  the  separatrix 
layer.  In  this  case,  there  would  also  be  a  hole  in  the  particle  scattering  plot  in  1  -  y  space 
which  is  reserved  for  regular  motion.  Fig.  10(b)  shows  a  instantaneous  plot  of  density  at 
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t  =  500ft-1.  It  is  interesting  to  notice  that  the  picture  is  quite  different  from  that  inferred 


above.  The  y-averaged  density  n  varies  linearly  in  x,  vanishes  at  the  position  about  one 
thermal  ion  gyroradius  from  the  wall,  as  seen  in  Theilhaber  and  Birdsall  simulation  plots  . 3 
This  observation  enables  us  to  surmise  that  the  motion  is  stochastic  in  the  whole  simulation 
region.  This  may  be  confirmed  in  the  scattering  plot  of  instantaneous  particle  position 
and  velocity  in  Figs.  10(c)-(d).  The  particles  near  the  center  of  the  vortex,  |x  -  2.5p,|  < 
0.5pj,  |y|  <  0.5pi,  axe  found  to  have  very  high  velocities  v  >  3.0t>7,.  The  motion  of  particles 
at  these  high  energies  (very  large  p,)  most  likely  are  stochastic  and  able  to  escape.  The 
increase  in  kinetic  energy  of  the  particles  near  the  bottom  of  the  vortex  well  is  a  consequence 
of  the  stochastic  motion,  allowing  individual  particle  to  cross  to  lower  potential  regions.  As 
the  total  particle  energy  is  conserved,  the  consequence  is  an  increase  in  kinetic  energy  and 
gyro-orbit  size  near  the  well  center  which  maintains  the  the  stochasticity. 

It  is  worth  pointing  out  that  during  the  simulation  we  observe  that  some  particles  axe 
trapped  near  the  source  region.  Thus,  some  fraction  of  the  particles  reintroduced  at  x  =  0 
as  source  particles  have  a  chance  to  be  trapped.  But,  compared  with  the  stochastic  flow,  the 
relative  number  of  trapped  to  escaping  particles  is  small.  By  comparing  particle  scattering 
plots  and  y-averaged  density  plots  at  t  =  250ft-1  and  t  =  1000ft-1 ,  we  confirm  the  approach 
to  quasi-steady-state  diffusive  flow. 
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Finally,  we  note  that  the  fully  developed  stochasticity  in  a  single  wave  found  in  Part 
I,  gives  an  estimate  for  the  diffusion  coefficient  of  D,  2:  0.07^.  This  is  sufficient  to  allow 
electrons  to  diffuse  across  the  sheath  region  at  a  rate  equal  to  that  of  the  fully  developed 
ion  diffusion  to  preserve  ambipolarity. 

9  Discussion  and  Conclusions 

Particle  simulation  of  motion  in  a  plasma  sheath  with  a  magnetic  field  perpendicular  to  the 
electric  field  indicates  there  is  quasi-steady-state  consisting  of  vortex  flow  on  the  scale  of 
the  sheath  dimension  and  a  small  amplitude  wave  spectrum  with  characteristic  wavelength 
of  the  order  of  the  ion  gyroradius.  There  is  a  continuous  transport  of  electrons  and  ions 
through  this  sheath,  even  in  the  absence  of  collisions,  from  the  plasma  source  to  the  sink 
at  the  material  wall. 

We  have  investigated  the  mechanism  for  the  transport  by  breaking  up  the  problem 
into  two  non- self- consistent  parts,  and  showing  that  they  can  combine  to  transport  both 
ions  and  electrons  across  the  sheath.  The  primary  driving  mechanism  is  the  resonance 
interaction  of  harmonics  of  the  vortex  frequency  with  the  ion  gyrofrequency  which  leads 
to  stochastic  diffusion.  For  the  parameters  of  the  vortex  size  and  strength,  and  the  self- 
consistent  ion  temperature,  the  numerical  integrations  of  single  particles  in  this  field  shows 
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that  the  stochastic  motion  transports  ions  across  the  sheath  at  a  rate  in  good  agreement 
with  that  found  in  the  simulation. 

Because  ti:  e  and  space  scales  are  not  well  separated,  analytic  techniques  using  secular 
perturbation  theory  are  not  sufficiently  accurate  to  quantitatively  calculate  the  stochas- 
ticity  and  resultant  diffusion.  They  do,  however,  uncover  the  basic  mechanisms  of  island 
formulation  and  destruction,  which  are  illustrated  by  surface  of  section  plots  of  the  single 
particle  motion.  Simple  analytic  estimates  of  the  diffusion  rate,  employing  the  concepts  of 
diffusion  theory  from  situations  in  which  the  scales  are  better  separated,  give  reasonable 
estimates  of  the  transport  rate. 

We  expect  that  the  stochastic  ion  motion  and  the  resulting  interception  of  ion  orbits  with 
walls  would  lead  to  macroscopic  charge  fluctuation  that  generate  a  wave  spectrum.  Indeed 
these  charge  fluctuations  and  waves  are  seen  in  the  simulation.  Without  inquiring  into 
the  details  of  the  spectrum,  we  show  that  even  a  single  wave,  of  the  frequency  and  wave 
number  corresponding  to  the  peak  observed  amplitude,  is  sufficient  to  resonate  with  the 
E  x  B  motion  of  the  electrons  in  the  vortex  to  generate  large  scale  stochasticity.  Additional 
waves,  at  the  same  overall  amplitude  level,  tend  to  make  the  stochasticity  more  uniform  with 
the  same  general  rate  of  transport.  For  E  x  B  motion  with  one  single  small  amplitude  wave 
there  is  better  quantitative  comparison  of  the  analytic  methods  with  the  single  particle 
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numerical  integrations.  This  is  because  the  wave  frequency  is  a  true  constant  of  motion 


which  made  the  system  well  suited  for  the  analytic  methods  used.  We  have  presented  these 
comparisons  and  shown  the  agreement  to  be  reasonable.  The  transport  rate  is  found  to  be 
sufficient  to  account  for  electron  loss  at  the  same  rate  as  the  ions,  and  thus  no  additional 
ambipolar  effects  need  be  postulated. 

The  parameters  explored  in  our  study  were  chosen  with  the  guidance  of  the  simulations, 
in  which  these  parameters  were  self-consistently  generated.  Significant  reduction  of  the  self- 
consistent  perturbation  parameters  or  of  the  gyroradius  led  to  considerably  less  stochasticity 
and  consequently  less  diffrs'  jn.  This  also  may  give  additional  insight  into  the  self-consistent 
problem,  as  it  in  „es  that  the  wave  strength  builds  up  until  ambipolarity  is  just  satisfied. 

The  fact  that  our  scaling  with  ^  (p,  and  vji)  in  Eq.  (60)  is  only  Bohm-like  for  fixed 
77,  *£•,  a-nd  £,  dose  not  mean  that  an  investigation  of  the  scaling  in  the  simulation  would 
not  give  Bohm-like  results;  i.e.,  as  parameters  T ,  and  B  change  in  the  simulation,  quantities 
a  and  b  (the  ratioes  of  the  vortex  size  to  the  gyro-radius)  could  either  remain  constant 
or  change  in  such  a  way  that  Eq.  (60)  would  give  the  Bohm-like  diffusion.  Our  work  points 
out  the  need  for  simulations  for  a  wide  parameter  range,  in  order  to  investigate  the  scaling, 
which  could  then  be  compared  to  the  type  of  calculation  we  have  performed. 

The  treatment  used  in  this  study  to  determine  the  transport  mechanism  and  calculate 
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the  rate  of  diffusive  loss  across  the  sheath  region  is  not  self  consistent.  In  fact  it  is  our 


ability  to  independently  vary  certain  parameters  and  investigate  the  consequence  of  the 
variation  that  gives  us  insight  into  the  transport  mechanism.  Nevertheless,  even  within  the 
conceptual  framework  of  our  study,  there  is  a  step  that  has  been  left  out,  the  determination 
of  the  waves  that  drive  the  electron  transport.  An  approach  to  this  problem  is  a  non-self- 
consistent  Fourier  analysis  of  the  ion  charge  fluctuations  within  the  context  of  ion  stochastic 
motion  and  transport.  However,  since  the  fluctuations  can  both  affect  the  ion  motion  and 
be  affected  by  the  electron  motion,  it  is  not  clear  that  the  resultant  fluctuation  spectrum 
will  closely  approximate  the  self-consistent  one.  In  any  case,  such  a  study  takes  us  beyond 
th^  framework  of  this  paper,  and  we  leave  it  for  future  work. 
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Figures 


Figure  1:  Characterization  of  the  physical  system  of  interest.  A  large  scale  coherent 
circular  flow  (or  vortex)  interacting  with  small  amplitude  traveling  waves. 

Figure  2:  Contours  of  the  potential,  Eq.  (4)  with  ip(u)  =  e  “7U,  with  5  =  1,7  =  2.5, 
a  =  1.25,  0  =  3.83,  v0o  =  0.5 vji,  Pi/L,  =  0,  and  kpi  =  3.83.  These  parameters  are 
chosen  to  compare  to  the  results  in  Refs.  2  and  3. 

Figure  3:  Plot  of  (x,y)  when  ui\t  =  2 irn,  n  =  1,2,3,...  (a)  (  =  0.015,  (b)  e  =  0.02. 

Figure  4:  Action  angle  coordinate  surface  of  section  plot  for  the  f  =  0.015  case.  Plot 
of  (J0,  60)  when  u qt  =  2nn,  n  =  1,2,3, ... 

Figure  5:  Surface  of  section  plot  with  four  waves.  Plot  of  ( x,y )  when  u>4t  =  2 xn, 
n  =  1,2,3,... 

Figure  6:  Four  test  particle  orbits  are  shown  just  at  the  interior  edge  of  the  vortex 
with  velocity  shear  S  =  0.015. 

Figure  7:(a)  Perspective  plot  of  the  potential  surface,  Eq.  (43)  with  jr-  =  jrp  = 
7[(  ~  iitrX  +  -  cos  $)J  7  =  ^  =  2.5,  a  =  1.59,  b  =  7.5  and  *  =  0.17. 

(b)  conto'  dot  of  0(x,  y)  in  (a). 
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Figure  8:  Plots  for  potential,  Eq.  (42).  The  vortex  center  is  located  at  x  =  0  and 


y  =  0.  (a)  Surface  of  section  for  N=12  particles,  ro  =  lp,,  voo  =  0.71UT-;  and  7  =  2.5; 
(b)  Single  particle  orbit  for  one  set  of  the  4-islands  in  (a). 

Figure  9:  Plots  for  potential,  Eq.  (43).  The  vortex  center  is  located  at  x  =  0  and 
y  =  0.  (a)  Surface  of  section  for  N  =  20  particles,  r0  =  0.5p,,  voo  =  1.55VT,,  and 
7  =  2.5;  (b)  Surface  of  section  for  N  =  20  particles, ro  =  0.5p,,u0o  =  1.7vr,-,  and 
7  =  2.5. 

Figure  10:  Typical  plots  for  diffusion  measurement.  The  vortex  center  is  located  at 
io  =  2.5p,,  y  =  0,  the  source  at  x  —  0  and  the  sink  at  x  =  5. Op,.  Here,  2500  particles 
are  advanced  in  advanced  in  time  by  500fi-1.  (a)  Number  of  particles  N,  escaping 
the  region  |x|  >  2xo  plotted  against  time  t.  Instantaneous  plots  at  t  —  500fi-1:  (b) 
Profile  of  the  y-averaged  ion  density;  (c)  ion  scatter  plot;  (d)  ion  velocity  plot. 
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